import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import rasterio as rio
from rasterio.plot import show
from rasterio.plot import show_hist
import os
import fiona

!pwd

/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis

!ls

Approx_SG.kml
debug.log
Displaying images.ipynb
Extracting_SG_from_Landsat8_images.ipynb
Generated DataFrames
HDB data
Landsat_Image_processing_final_experiments.ipynb
Landsat8_image_processing_initial_experimentation.ipynb
Population_data_World_Bank.ipynb
Population_density_and_Housing_density_analysis.ipynb
Predicting the future of SG.ipynb
SG_cropped_Landsat8_images
SG_Landsat8_images
SG_shapefiles

os.chdir('SG_cropped_Landsat8_images/')

!pwd

/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis/SG_cropped_Landsat8_images

!ls

2013
2014
2015
2016
2017
2018
2019
2020

os.chdir('2013')

SG_Bands_2013 = {}

for i in range(1,12):
    var = 'B'+str(i)
    
    img_file = 'SGimgB' + str(i) + '.TIF'
    SG_Bands_2013[var] = rio.open(img_file, 'r')

os.chdir('..')

os.chdir('2014')

SG_Bands_2014 = {}

for i in range(1,12):
    var = 'B'+str(i)
    
    img_file = 'SGimgB' + str(i) + '.TIF'
    SG_Bands_2014[var] = rio.open(img_file, 'r')

os.chdir('..')

os.chdir('2015')

SG_Bands_2015 = {}

for i in range(1,12):
    var = 'B'+str(i)
    
    img_file = 'SGimgB' + str(i) + '.TIF'
    SG_Bands_2015[var] = rio.open(img_file, 'r')

os.chdir('..')

os.chdir('2016')

SG_Bands_2016 = {}

for i in range(1,12):
    var = 'B'+str(i)
    
    img_file = 'SGimgB' + str(i) + '.TIF'
    SG_Bands_2016[var] = rio.open(img_file, 'r')

os.chdir('..')

os.chdir('2017')

SG_Bands_2017 = {}

for i in range(1,12):
    var = 'B'+str(i)
    
    img_file = 'SGimgB' + str(i) + '.TIF'
    SG_Bands_2017[var] = rio.open(img_file, 'r')

os.chdir('..')

os.chdir('2018')

SG_Bands_2018 = {}

for i in range(1,12):
    var = 'B'+str(i)
    
    img_file = 'SGimgB' + str(i) + '.TIF'
    SG_Bands_2018[var] = rio.open(img_file, 'r')

os.chdir('..')

os.chdir('2019')

SG_Bands_2019 = {}

for i in range(1,12):
    var = 'B'+str(i)
    
    img_file = 'SGimgB' + str(i) + '.TIF'
    SG_Bands_2019[var] = rio.open(img_file, 'r')

os.chdir('..')

os.chdir('2020')

SG_Bands_2020 = {}

for i in range(1,12):
    var = 'B'+str(i)
        
    img_file = 'SGimgB' + str(i) + '.TIF'
    SG_Bands_2020[var] = rio.open(img_file, 'r')

os.chdir('..')

Function to display Bands 1 - 9 for each of the cropped Landsat Image

def show_bands(dictionary):
    
    Bands = dictionary
    fig, ax = plt.subplots(3,3, figsize = (15,10))

    num = 0
    for i in range(3):
        for j in range(3):
            num += 1
            key = 'B' + str(num)
        
            ax[i,j].imshow(Bands[key].read(1))
            ax[i,j].set_xlabel(key)
        

Experimental function to check which bands would return best difference between Land and Water

• Experimenting to see which combination of Bands produces the best differentiation between Land and Water while avoiding clouds.
• This combination may be different for each Landsat scene.
• Using the cropped image of SG allows for much faster experimentation.

def experiment(dictionary, band_num):
    
    Bands = dictionary
    main_key = 'B' + str(band_num)
    
    fig, ax = plt.subplots(3,3, figsize = (15,10))
    
    num = 0
    for i in range(3):
        for j in range(3):
            num += 1
            key = 'B' + str(num)
            
            if num == 8: # Excluding Band 8 since it has higher resolution and cannot be subtracted due to shape mismatch
                continue
                
            ax[i,j].imshow(Bands[main_key].read(1) - Bands[key].read(1))
            ax[i,j].set_xlabel(main_key + " - " + key)
    

show_bands(SG_Bands_2013)

experiment(SG_Bands_2013, 5)

plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2013['B5'].read(1) - SG_Bands_2013['B4'].read(1) )

<matplotlib.image.AxesImage at 0x1f7221848b0>

show_bands(SG_Bands_2014)

experiment(SG_Bands_2014, 5)

plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2014['B5'].read(1) - SG_Bands_2014['B3'].read(1) )

<matplotlib.image.AxesImage at 0x1f7210e8160>

show_bands(SG_Bands_2015)

experiment(SG_Bands_2015, 5)

plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2015['B5'].read(1) - SG_Bands_2015['B3'].read(1) )

<matplotlib.image.AxesImage at 0x1f7221d2e50>

show_bands(SG_Bands_2016)

experiment(SG_Bands_2016, 5)

plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2016['B5'].read(1) - SG_Bands_2016['B3'].read(1) )

<matplotlib.image.AxesImage at 0x1f71f821d90>

show_bands(SG_Bands_2017)

experiment(SG_Bands_2017, 5)

plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2017['B5'].read(1) - SG_Bands_2017['B3'].read(1) )

<matplotlib.image.AxesImage at 0x1f71f70c490>

show_bands(SG_Bands_2018)

experiment(SG_Bands_2018, 5)

plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2018['B5'].read(1) - SG_Bands_2018['B2'].read(1) )

<matplotlib.image.AxesImage at 0x1f71d45d490>

show_bands(SG_Bands_2019)

experiment(SG_Bands_2019, 5)

plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2019['B5'].read(1) - SG_Bands_2019['B3'].read(1) )

<matplotlib.image.AxesImage at 0x1f71cddf0d0>

show_bands(SG_Bands_2020)

experiment(SG_Bands_2020, 5)

plt.figure(figsize = (15,10))
plt.imshow(SG_Bands_2020['B5'].read(1) - SG_Bands_2020['B3'].read(1) )

<matplotlib.image.AxesImage at 0x1f71f587910>

SG_imgs = {}
SG_imgs[2013] = SG_Bands_2013['B5'].read(1) - SG_Bands_2013['B3'].read(1)
SG_imgs[2014] = SG_Bands_2014['B5'].read(1) - SG_Bands_2014['B3'].read(1)
SG_imgs[2015] = SG_Bands_2015['B5'].read(1) - SG_Bands_2015['B3'].read(1)
SG_imgs[2016] = SG_Bands_2016['B5'].read(1) - SG_Bands_2016['B3'].read(1)
SG_imgs[2017] = SG_Bands_2017['B5'].read(1) - SG_Bands_2017['B3'].read(1)
SG_imgs[2018] = SG_Bands_2018['B5'].read(1) - SG_Bands_2018['B2'].read(1)
SG_imgs[2019] = SG_Bands_2019['B5'].read(1) - SG_Bands_2019['B3'].read(1)
SG_imgs[2020] = SG_Bands_2020['B5'].read(1) - SG_Bands_2020['B3'].read(1)

# Images that differentiate land and water the best

fig, ax = plt.subplots(len(SG_imgs)//2,2, figsize = (15,15))

num = 2013

for i in range(len(SG_imgs)//2):
    for j in range(2):
        ax[i, j].imshow(SG_imgs[num])
        ax[i, j].set_xlabel(num)
        num += 1
        
plt.tight_layout()

# Images flattened to 1D arrays
SG_imgs_flattened = {}
for i in range(2013,2021):
    print(SG_imgs[i].shape)
    SG_imgs_flattened[i] = SG_imgs[i].flatten()
    print(SG_imgs_flattened[i])

(1216, 1877)
[0 0 0 ... 0 0 0]
(1216, 1877)
[0 0 0 ... 0 0 0]
(1216, 1877)
[0 0 0 ... 0 0 0]
(1216, 1877)
[0 0 0 ... 0 0 0]
(1216, 1877)
[0 0 0 ... 0 0 0]
(1216, 1877)
[0 0 0 ... 0 0 0]
(1216, 1877)
[0 0 0 ... 0 0 0]
(1216, 1877)
[0 0 0 ... 0 0 0]

Plotting histograms to visualize land and water demarcation

• Most of the land and water pixels are clearly segregated.

fig, ax = plt.subplots(len(SG_imgs)//2,2, figsize = (15,15))

num = 2013

for i in range(len(SG_imgs)//2):
    for j in range(2):
        ax[i, j].hist(SG_imgs_flattened[num], bins = 100, color = 'r')
        ax[i, j].set_xlabel("Histogram for " + str(num))
        num += 1
        
plt.tight_layout()

Cropped Landsat images of SG with respective Histograms

fig, ax = plt.subplots(len(SG_imgs),2, figsize = (15,15))

num = 2013

for i in range(len(SG_imgs)):
    ax[i,0].imshow(SG_imgs[num])
    ax[i,1].hist(SG_imgs_flattened[num], bins = 100, color = 'r')
    ax[i,0].set_xlabel(num)
    ax[i,1].set_xlabel("Histogram for " + str(num))
    num += 1
        
plt.tight_layout()

#show_histogram(SG_Bands_2013)

# Finding the total number of pixels in the manually drawn shapefile

SG_Bands_2013_flat = SG_Bands_2013['B9'].read(1).flatten()
#type(SG_Bands_2013_flat)

total_pixels_in_cropped_img = (SG_Bands_2013_flat > 0).sum()
print(total_pixels_in_cropped_img)

1372139

# Thresholding
def Thresholding():
    SG_binary_imgs_flat = {}
    for i in range(2013, 2021):
    
        SG_binary_imgs_flat[i] = np.zeros(len(SG_imgs_flattened[i]))
        
        for j in range(len(SG_imgs_flattened[i])):
            if SG_imgs_flattened[i][j] > 20000:
                SG_binary_imgs_flat[i][j] = 1
    
    return SG_binary_imgs_flat
                

SG_binary_imgs_flat = Thresholding()

SG_binary_imgs_flat_copy = SG_binary_imgs_flat.copy()

# Creating Binary image

SG_binary_imgs = {}
for i in range(2013, 2021):
    SG_binary_imgs[i] = SG_binary_imgs_flat_copy[i].reshape((1216,1877))

Binary Images

fig, ax = plt.subplots(len(SG_imgs)//2,2, figsize = (15,15))

num = 2013

for i in range(len(SG_imgs)//2):
    for j in range(2):
        ax[i, j].imshow(SG_binary_imgs[num])
        ax[i, j].set_xlabel("Binary Image: " + str(num))
        num += 1
        
plt.tight_layout()

#Counting number of pixels of water

SG_water_pixels = {}
for i in range(2013, 2021):
    SG_water_pixels[i] = (SG_binary_imgs_flat[i] == 1).sum()
    

#Counting number of pixels of land

SG_land_pixels = {}
for i in range(2013, 2021):
    SG_land_pixels[i] = total_pixels_in_cropped_img - SG_water_pixels[i]

for i in range(2013, 2021):
    print("Year: " + str(i))
    print("     #Pixels: " + str(SG_land_pixels[i]))

Year: 2013
     #Pixels: 774331
Year: 2014
     #Pixels: 800777
Year: 2015
     #Pixels: 791720
Year: 2016
     #Pixels: 795768
Year: 2017
     #Pixels: 801289
Year: 2018
     #Pixels: 790573
Year: 2019
     #Pixels: 798514
Year: 2020
     #Pixels: 798597

# Dictionary that has the area of SG Land in squared km based on no. land pixels
SG_area = {}

 # Data from https://data.gov.sg/

Formula converting number of pixels to land area

$C = constant$

$C * NumPixels = Area$

From 2015 data (Data from https://data.gov.sg/),

$C * 791720 = 719.1$

$ -> C = 0.0008947559682658218 $

C = 719.1/791720

for i in range(2013, 2021):
    SG_area[i] = [C * SG_land_pixels[i]]

719.1/791720

0.0009082756530086394

print(SG_area)

{2013: [703.3059946698328], 2014: [727.3262525892992], 2015: [719.1], 2016: [722.776699843379], 2017: [727.7912897236397], 2018: [718.058207825999], 2019: [725.2708247865406], 2020: [725.3462116657404]}

!pwd

/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis/SG_cropped_Landsat8_images

os.chdir('..')

os.chdir('HDB data')

!ls

dwelling-units-under-hdbs-management-by-town-and-flat-type.csv
estimated-percentage-of-singapore-resident-population-in-hdb-flats.csv
estimated-resident-population-in-hdb-flats-by-town.csv
Houses_by_type_2020.png
metadata-estimated-resident-population-living-in-hdb-flats.txt
metadata-number-of-residential-units-under-hdb-s-management.txt
metadata-total-land-area-of-singapore.txt
SG_population_data.csv
total-land-area-of-singapore.csv
World_bank_population_data.csv

SG_land_area_data = pd.read_csv('total-land-area-of-singapore.csv')

SG_land_area_data

	year	total_land_area
0	1960	581.5
1	1961	581.5
2	1962	581.5
3	1963	581.5
4	1964	581.5
5	1965	581.5
6	1966	581.5
7	1967	583.0
8	1968	584.3
9	1969	585.3
10	1970	586.4
11	1971	586.4
12	1972	586.4
13	1973	586.4
14	1974	587.6
15	1975	596.8
16	1976	602.0
17	1977	616.3
18	1978	616.3
19	1979	617.8
20	1980	617.8
21	1981	617.9
22	1982	618.1
23	1983	618.1
24	1984	620.2
25	1985	620.5
26	1986	621.7
27	1987	622.6
28	1988	625.6
29	1989	626.4
30	1990	633.0
31	1991	639.1
32	1992	641.0
33	1993	641.4
34	1994	646.1
35	1995	647.5
36	1996	647.5
37	1997	647.8
38	1998	648.1
39	1999	659.9
40	2000	682.7
41	2001	682.3
42	2002	687.1
43	2003	693.4
44	2004	696.2
45	2005	697.9
46	2006	699.5
47	2007	705.1
48	2008	710.2
49	2009	710.3
50	2010	712.4
51	2011	714.3
52	2012	715.8

SG_land_area = pd.DataFrame(SG_area, index = ["total_land_area"])

SG_land_area

	2013	2014	2015	2016	2017	2018	2019	2020
total_land_area	703.305995	727.326253	719.1	722.7767	727.79129	718.058208	725.270825	725.346212

SG_land_area = SG_land_area.transpose()

SG_land_area = SG_land_area.reset_index()

SG_land_area = SG_land_area.rename(columns ={"index": "year"} )

Dataframe that tabulates total land area of SG over the lears

SG_land_area_data = SG_land_area_data.append(SG_land_area)

SG_land_area_data

	year	total_land_area
0	1960	581.500000
1	1961	581.500000
2	1962	581.500000
3	1963	581.500000
4	1964	581.500000
...	...	...
3	2016	722.776700
4	2017	727.791290
5	2018	718.058208
6	2019	725.270825
7	2020	725.346212

61 rows × 2 columns

plt.rcParams.update({'font.size': 22})
plt.figure(figsize = (20,7))
plt.plot(SG_land_area_data['year'], SG_land_area_data['total_land_area'], '-o', label = 'Actual Area')
plt.plot(SG_land_area['year'], SG_land_area['total_land_area'], '-or', label = 'Estimated Area from Landsat 8 images')
plt.xlabel("Year")
plt.ylabel("Area in square km")
plt.legend()
plt.grid()

os.chdir('..')

os.chdir('Generated DataFrames/')

!pwd

/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis/Generated DataFrames

SG_land_area_data.to_csv('SG_land_area_data.csv', index= False)